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Abstract 

Computer algebra procedures to manipulate pseudo-differential operators are im- 
plemented to perform calculations with integrable models. We use lazy evaluation 
and streams to represent and operate with pseudo-differential operators. No order 
of truncation is needed since terms are produced on demand. We give a series of 
concrete examples using the computer algebra language MAPLE. 



*Email: brunclli@fsc.ufsc.br 



1 



1 Introduction 



In many calculations in physics we are faced with the problem of obtaining the product of 
infinite power series such as 

oo oo 

This is usually done through the convolution formula 

oo 3=i 

F^G = 5^x^5^/,^.„,. (2) 

i=0 j=0 

Implementing these products of series in a computer algebra system is usually done keep- 
ing track of indexes and expansion order of truncation, which is very cumbersome. The 
algebraic language MAPLE provides the package powseries which allows us to manipu- 
late formal power series without the need to specify the truncation order. We can ask as 
many terms as desired. Encapsulated in this package is a powerful paradigm of modern 
computation, streams with delayed or lazy evaluation, or simply streams for short. Streams 
are natural data structures to implement programs which use infinite objects such as se- 
quences, series, etc. P and, as we will see, pseudo-differential operators (PDO). The use 
of streams frees us of the task of taking care of expasion order of truncation and avoids the 
proliferation of summation indexes. 

PDO are widely used in the theory of integrable models [SElll]- Calculations with 
them are systematic, such as obtaining equations of motion via Lax pairs as well the gen- 
eration of conserved charges via fractional powers of the Lax operators. Therefore, these 
calculations can be performed by the assistance of some sort of computer algebra system. 
In fact, a powerful program^ aimed to this task was introduced in [6J using the language 
AXIOM |Zj (formerly known as Scratchpad). AXIOM was used because stream is an prim- 
itive data structure in this language and the existence of a large library of "domains" allow 
us to construct very generic algorithms. Unfortunately, AXIOM is not a widespread lan- 
guage and is not commercially available any more (efforts are being made to make AXIOM 
a public domain software, see http://www.nongnu.org/axiom/). In this paper we will 
provide such algorithm in MAPLE, however, we will have to construct streams with lazy 
evaluation form scratch. This is in fact very straightforward jHIH]^ and in doing so we hope 
to make these beautiful and elegant techniques more popular to physicists (see |9| for a nice 
apology on these techniques). 

At this point we would like to make some general comments. Laziness evaluates pro- 
cedure arguments only when the procedure uses them, otherwise the evaluation is delayed. 
This make possible to manipulate infinite data structures, such as the series and PDOs 
in a simple way. This technique of evaluation is the key ingredient of our algorithms. We 
will also use a functional programming style as much as possible to construct our algorithms 

^See also |S] for other programs dealing with integrable models. 

^Let us point out that the paper |H| was written for a version of MAPLE where nested lexical scoping 
was not available and therefore it had to be simulated. 
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in MAPLE. This will make our programs small and almost free of loops, special cases state- 
ments, synchronization of expansion orders, assignments, etc. In this style of programming, 
mathematical statements like 

0! = 1, 

n\ = n{n — 1)! 

could be written, for instance in MAPLE, simply as 
> new_f act :=n-> ' if ' (n=0 , 1 ,n*new_f act (n-1) ) ; 

an almost literal transcription of a mathematical "formula" . This should be contrasted with 
a similar program written in an imperative programming style. In functional programming 
the code of programs is compact and elegant, easy to read and modify. In this example 
new_f act is a recursive function or procedure and they have a key role in the functional style 
of programming. Also, new_f act has a terminal clause, n = 0, and we will see that thanks 
to laziness it will be possible to construct recursive functions without terminal clauses, which 
we call corecursive functions. In the functional style of programming procedures must be 
first class, also called higher order procedures. That means that procedures can accept other 
procedures as parameters or can produce another procedure as their result. Laziness and 
the functional style of programming are some of the main characteristics of modern (pure) 
functional programming languages. Is the simultaneous use of them that makes possible 
to describe abstract mathematical objects such as PDO in a concrete way. Despite the 
fact that MAPLE is not a functional language, such as Haskell jTHj, procedures are first 
class and we can implement laziness evaluation and streams while taking advantage of the 
resources available in this powerful and widely used language. Even though programs using 
lazy evaluation (or programs written in functional languages in general), from a semantic 
point of view, are compact and have elegant code most of these programs use recursive calls 
or recursive defining relations that causes recomputation of terms several times. In order 
to minimize this drawback we will use the MAPLE remember option in our procedures^, in 
that way a procedure stores its results in a remember table. Of course, this costs memory 
and beyond lower orders we can get memory overflow. 

This paper is organized as follows. In sections 2 and 3 we introduce with detail the 
general procedures to implement streams with lazy evaluation. We hope these can motivate 
readers to use these techniques in their own programs. In section 4 we review facts about 
PDO to be used in section 5 about integrable models. In section 6 we present our program 
and its routines. A serie of concrete examples is given in section 7. We end the paper with 
some conclusions in section 8. 

2 Delay and Force 

We implement delayed evaluation following [T] closely. The procedure^ delay packs any ex- 
pression in a way that it can be evaluated later on demand. We use a MAPLE module to rep- 

•^For self-referencial procedures we can also use a fixed point method |H1 ■ 

■^Also, following P we adopt descriptive names for procedures, parameters, etc, trying to make programs 
as self-documenting as possible. 



resent a record with two exports: delayed. object and forced?. In the delayed. object 

export we use the body of a procedure of no arguments as the delayed expression. The 
procedure force will execute the procedure, forcing the evaluation. Also, this vahie will 
be stored and subsequent forcing of the same delayed object will return the stored value 
avoiding inefficiency problems in recursive programs with streams. The export forced? 
keeps track if the expression was already forced. 

> delay : =proc (s : :uneval) 

option remember; 
local delay_record; 
delay_record : =module () 

export delayed_object .forced? ; 
end module; 

delay_record: -forced? : =f alse ; 
delay_record : -delayed_ob j ect : =proc () 

eval(s) ; 
end proc; 

Delay (delay_record) ; 
end proc: 

> f orce:=proc(s) 

option remember; 

local delay_record; 

delay_record : =op ( 1 , s) ; 

if not delay_record: -forced? then 

delay_record: -delayed_object : =delay_record: -delayed_object () ; 

delay_record : -forced? : =true ; 
end if ; 

delay_record : -delayed_ob j ect ; 
end proc: 

> 'print /Delay' := procO 

"Delayed" ; 
end proc: 

As an simple example consider 

> de:=delay(diff (x"2,x)) ; 

de := "Delayed" 

> eval(op( [1,1,2] ,de)) ; 

false 
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> f orce(de) ; 



2x 

> eval(op( [1,1,2] ,de)) ; 

true 

Usually, the arguments of a procedure are evaluated before the procedure is invoked. 
This is known as applicative-order evaluation. This behavior will return an error in the 
following procedure 

> ab:=proc(a,b) 

'if '(a=0,l,b); 
end proc: 

> ab(0,l/sin(0)) ; 

Error, numeric exception: division by zero 

Using delay we can postpone the evaluation of the argument of a procedure until it is 
needed. This is called normal-order evaluation or lazy evaluation. In this way 

> ab(0,delay(l/sin(0))) ; 

1 

> ab(l,delay(l/sin(0))) ; 

"Delayed" 

and no error is produced as long as we do not force the evaluation. This example shows 
us that we can perform useful computations even when values of some parameters would 
produce errors or are not known. We will exploit lazy evaluation in the next section to 
construct infinite structures. 

3 Streams 

Stream is a data structure representing a sequence of infinite terms. We implement streams 
as lists with delayed evaluation. Following the usual LISP terminology we introduce the 
constructor cons_streaiii and two selectors, car_streaiii and cdr_streaiii J^. The proce- 
dure cons_stream takes the two arguments head and tail and returns a compound object, 
a delayed list with the arguments as its parts 



> cons_stream:=proc(head: :uneval, tail : :uneval) 

delay ( [head, tail] ) ; 
end proc: 

With cons_stream we can form pairs whose elements are pairs, and so on, i.e. we can 
construct infinite Usts. Given a stream we can extract its head or first element using the 
selector car_streain 

> car_stream: =proc (s) 

op(l,f orce(s)) ; 
end proc: 

The selector cdr_stream selects the sublist consisting of all but the first item 

> cdr_streain:=proc(s) 

op(2,force(s)) ; 
end proc: 

Noticing that nested applications of car_stream and cdr_stream} allows us to extract any 
element of a stream we introduce the procedure ref _stream which takes as arguments a 
stream s and a number n and returns the item of the stream. It is customary to number 
the elements of the stream beginning with 0. The method for computing ref .stream is 
the following: for n — 0, ref .stream should return the head of the stream, otherwise, it 
should return the {n — 1) item of the tail of the stream. Therefore, we get the recursive 
procedure 

> ref .stream :=proc(s,n) 

if n=0 then 

car_stream(s) ; 

else 

ref _stream(cdr_stream(s) ,n-l) ; 
end if ; 
end proc: 

Now we can construct our first infinite structure, a infinite list of ones 

> stream_of _ones : =cons_stream(l , stream_of _ones) ; 

strearri-of -Ones := "Delayed" 
The stream stream_of _ones creates a list whose head equals 1 

> car_stream(stream_of _ones) ; 

1 

but whose tail is not computed 



> cdr_stream(streain_of _ones) ; 

"Delayed" 

> car_stream(cdr_stream(stream_of _ones) ) ; 

1 

> ref _streain(stream_of _ones , 100) ; 

1 

Note that streams_of _ones is corecursive and thanks to the laziness the recursive evalu- 
ation is avoided. The main point to be notice with this example is that due to laziness a 
recursive equation has become an algorithm. 

Given a function a{k) that allow us to compute the k^^ element of a stream we can 
generate the stream, starting from the n*^ element, a{n + 0), a{n + 1), a{n + 2), . . . by the 
procedure make_stream 

> make_stream:=proc(a,n) 

option remember; 

cons_stream(a(n) ,make_stream(a,n+l) ) : 
end proc: 

In order to display the first m terms of a stream s (or the m — 1 first items ) we use the 
procedure display_streain 

> display_stream: =proc (s ,m) 

if m=0 then 
NULL; 

else 

car_streain(s) ,display_streain(cdr_streain(s) ,m-l) ; 
end if ; 
end proc: 

Now we can represent the series (pQ) as streams 

> F:=make_stream(proc(n) 

local f; 

printf ("V/.a.",n) ; 
'if ' (n=0,f [0] ,f [n]*x"n) ; 
end proCjO) : 

> G:=make_stream(proc(n) 

local g; 

printf ("V/.a.",n) ; 
'if '(n=0,g[0] ,g[n]*x-n); 
end proCjO) : 
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Note that when a term is computed 

> display_stream(F, 10) ; 

0.1.2.3.4.5.6.7.8.9. 

/o, fix, f2X^, /3X^ /4X^, /5X^ /e /t x'^, /s X^ /gX^ 

it is never recomputed 

> display_streaiii(F, 11) ; 

10. 

/o, fix, f2X^, fiX^, /4X^, /5X^ /eX^ /7X^ /s X^ /gX^ /lo X^° 

Let us represent the infinite series in as streams in the form of the head /o com- 
pounded with the tail power series F by cons_stream, which we represent by " : " 

F = fo:F. 

The sum of the series F and G can then be represented by 

F(BG={fo:F)(B{go:G) = (/o + 9o) : (F ®G) . 
The product of a series F by a constant c is simply 

c-F = c-(/o:F) = (c/o):(c-F). 
Finally, the product of the series F and G, as in can be written as 

F^G=ifo:F)^ig,:G) = ifogo): [{fo -G) ® (F ^ G)] . 
We can code in MAPLE all these algorithms easily as 

> add_stream:=proc(sl,s2) 

cons_stream(car_streain(sl)+car_streain(s2) , 

add_streaiii(cdr_streain(sl) , cdr_streani(s2) ) 

); 

end proc: 

> display_streaiii(add_streaiii(F,G) , 10) ; 
0.1.2.3.4.5.6.7.8.9. 

/o + 5-0, fix + gi X, /2 x^ + g2 x^ /s x^ + 5(3 x^ f^ x^ + g^ x^, /s x^ + g^ x^ 
/e x^ + gQ x^ /7 x'^ + gj x\ /§ x^ + 5-8 /g x^ + 5(9 x^ 
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> scalar_multiply_streaiii: =proc (constant , s) 

cons_stream( constant *car_stream(s) , 

scalar_multiply_streain(constant , cdr_streaiii(s) ) 

); 

end proc: 

> display_stream(scalar_multiply_stream(c,F) , 10) ; 



c/o, cfix, C/2X^, cfsx^, cfiX^, cf^x^, cfex^, cfjx'^, cfsx^, cfgx^ 

> multiply_stream: =proc(sl , s2) 

cons_stream(car_streain(sl) *car_stream(s2) , 
add_streain( 

scalar_multiply_streain(car_stream(sl) , cdr_streaiii(s2) ) , 
multiply_streain(cdr_stream(sl) , s2) 

) 

); 

end proc: 

> display_streani(multiply_stream(F .G) ,5) ; 

/o 5-0, fo9ix + fi go x, fo g2 x'^ + figix + /2 5-0 a;^ /o ga + /i 5-2 + /2 gi x^ + fs go 

/o g4 x'^ + /i 5-3 + /2 5-2 x'^ + /a 5(1 + f^go x^ 

Now we see how the implementation of Q using streams is free of indexing and order 
of truncation "bureaucracy". We can go on and define the quotient, reversion and other 
operations on power series (see references jHl 1121 CSl Cll IIS]^)- 

A very useful operation is to apply a transformation to each element in a stream and 
generate the stream of results. The procedure map_stream takes as arguments a procedure 
/ of one argument and a stream s, and returns a stream of the results produced by applying 
the procedure to each element in the stream 

> map_stream: =proc (f , s) 

cons_stream(f (car_stream(s) ) ,map_stream(f , cdr_stream(s) ) ) ; 
end proc: 

> display_streani(map_stream(term->dif f (term,x) ,F) ,10) ; 



0, /i, 2/2X, 3/3X2, 4/4x3, 5/5X^ 6/6X^ 7^X^ S/gX^ g/gX^ 

Similar to map_stream we define zip_stream which takes n streams and maps a function 
/ with n arguments onto them 



^See especially ^21 ) where in a pure functional language, such as Haskell, the "music of streams is 
played" or ^3] where the " functional pearls shine" . 
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> zip_stream: =proc (f , streams) 

cons_stream(f (map(car_stream, [args [2 . . -1] ])[]), 
zip_stream(f ,map(cdr_stream, [args [2 . . -1] ] ) [] ) 

); 

end proc: 



> display_streani(zip_stream((x,y)->x"2+y"2,F,G) ,5) ; 



fo' + go', fi' + g,' x\ + g,' x\ x' + gs' x^ f,' x' + g,' x' 



The procedure f ilter_streain selects some elements of a given stream s forming a new 
stream which only contains those elements for which the given predicate predicate? is 
true 

> f ilter_stream : =proc (predicate? , s) 

if predicate? (car_streain(s)) then 



else 

f ilter_streain(predicate? , (cdr_stream(s) ) ) ; 
end if ; 
end proc: 

> display_streaiii(f ilter_stream(term->isprime(degree(term,x)) ,F) , 10) ; 



Using procedures such as map_stream, zip_stream and f ilter_stream we can im- 
plement algorithms in a sort of signal-flow structure ( "signals" that flow from one "process 
stage" to the next) if we use streams to represent "signals" and stream operations to im- 
plement each of the "process stages" [T]. 

4 Pseudo Differential Operators 

A differential operator is the finite sum 



cons_stream(car_streain(s) , 

filter_streaiii (predicate? , cdr_stream(s) ) 

); 



11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28.29.30.31. 



f2x\ /3X^ /5X^ /iix", /l3X^^ /l7a;^^ fi9X^\ /23a;^^ /29a;' 



,29 



n 




(3) 



i=0 



where the coefficients Pi[u] are differential functions and 




in 



i.e., d is the differential operator with respect to the variable x. Here we are considering a 
ring of differential operators in one independent variable x and one dependent variable u. 
For any ^, j > we have 

and for a differential function Q we have the Leibniz rule 

fc=0 ^ ^ 

where 

g(fc) ^ ^ _ 

These two rules allow us to define the product of any two differential operators making the 
space of all differential operators a noncommutative ring with the constant function 1 as 
the identity multiplication operator. 

Since (jS)) is a polynomial in the total derivative d with differential functions as co- 
efficients we obtain pseudo-differential operators (analogous of Laurent series) if we allow 
negative powers of d. In this way a formal pseudo-differential operator (PDO) is the infinite 
series 

n 

P= Yl (6) 

L has order n if its leading coefficient is not identically zero and by convention, the zero 
PDO is said to have order — oo. 

For PDOs the relation ^ is still valid, however, the Leibniz rule (0) generalizes to 

guQ ^ J2 n{n-l)---in-k + l) g(,)^„-fc ^ 



k=0 



which is valid to any n and reduces to (0) if > 0. For n > ^ can be written as 

d-^Q = ^(-1)' ( ^ ^ ~ M . 



Again, Q and (|7j) allow us to define the product of any two PDOs. However, for algorithm 
implementation we will use the following result JH] (see also |17j, exercise 5.20). For the 
PDOs 



P = P[u,d]= P^M^\ 

i=— oo 
m 

Q = Q[u,d] = QjM9' : 



j=-oo 
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we can associate the formal Laurent series 

n 

P[u,p\ = J2 P^MP'^ 
i=—oo 
m 

j = -OD 

by formally substituting the operator d for the variable p. The product of P and Q, 

R[u, d]=P-Q , 
is then determined by the Laurent series 

Id'Pd'Q 



1=0 ^ 



Rewriting Q in the form 



p = J2p,[u,p], 



1=0 

oo 



Q = YQj[u,p] 

j=0 



the expression (jiup assumes the form 

oo k k—i 



1 d'Pj d'Q 

fc=0 4=0 j=0 

As an example consider 



fc=0 4=0 j=0 ^ 



then 



p = d + u^d^^ 
Q = d^ + ud^' 



P[u,p] = p + v^p ^ , 
Q[u,p] = p^ + up~^ , 



and (fTTj) yields 



R[u,p] = PoQo + PoQi + PiQo + P0Q2 + PiQi + P2Q0 



dp dx dp dx dp^ dx^ 

+ -U + V^p + V^p"'^ + UxP'^ — U^UxP~^ + U^UxxP~'^ + • • • 
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and therefore 



It can be shown that every nonzero PDO P has an inverse and every PDO of order 
n > 0, with coefficient P^, has an n*^ root 

N=^/P= + ao + a_id-^ + a^2d'^ + a^g^-^ + • • • . (12) 



Imposing 

{NY = P , 

we obtain a linear system of equations for the coefficients a„ of which can be recursively 
solved. Also, fractional powers 

pWn = (^)- (13) 
can be defined. As an example, for the operator 

L = d'^ + u, (14) 

we have 

3 3 1 1 
L^/-^ = d^ + -ud+-u^ + -{3u^ + u^x)d-^ -—{6uu^ + u^^^)d-^ + -- - . (15) 

5 Integrable Models 

Integrable systems have been studied extensively from various points of view j2l El E] after 
the discovery of the inverse scattering transform. Now it is well known that integrable 
systems of the form 

ut = K[u] , 

where u = u{x,t) and i^' is a differential function of u, can be represented in the form of 
the so called Lax equation [iHl CEl EH] 

^ = [5,^], (16) 

where the Lax pair L and B are in general PDO with scalar coefficient functions. For 
instance, the generalized KdV hierarchy is described in terms of the Lax pair of order n 

1201 ini 1221 23J 



■2 , 



B={L'J-)_,, k^n, (17) 

where Uo{x,t), ■ ■ ■ ,Un-3{x,t),Un-2{x,t) represent the dynamical variables of the system 
and ( )^ denotes the part of the PDO with nonnegative powers. In more general cases 
(nonstandard Lax representation) we can have 

B = {)>,, (18) 



representing the part of a PDO with terms 9", n > d. In this way (^L^^^^ and (^L^.^^ 
are equivalent. The evolution equations, or the k^^ flow of the hierarchy, are given by 

9L„ 



>o 



dtk 



[{L'Jn^^L^]- (19) 



The scalar Lax equation ()16|) is the compatibility condition for the system of linear 
equations 

Lip = X'lp , 

imposing the constancy in time of the spectral parameter A. 

The KdV hierarchy is obtained from ()17|1 setting n = 2 and the KdV equation corre- 
sponds to = 3 (with «_! = 0, M = Mo). Therefore (see (ITH|) ) 

L = (9^ + M , 

B={L'/^)^ = d'' + '^{du + ud) , (20) 
and the Lax equation (|TH|l or |T!I|) yields the KdV equation 

Ut = -^U^xx + -jUU^ . (21) 

Given a PDO P = Yli=-oo ^i'^i residue is deflned as the coefficient of d^^, 

ResP = P_i, (22) 

and its trace is deflned as 

TrP= / dxResP= / dx P_i . 



It is easy to show that Tr satisfles the cyclicity property |2T 

Tr Pg = TrQP , 

which implies 



Tr[P,g] = 0, (23) 



for any two PDO P and Q. 

From the Lax equation ()19|) we can write 



^ = [W")+,^n/l, (24) 



for arbitrary m. Now, it can be shown using (|19p and (j24|) that distinct flows commute 

[221123] 

dtmdtk dtkdtm 
Also, taking the trace of and using we get 

|-Tr (L-/") = . 

Thus, we define the conserved charges, under any flow, as 

H^= I dxhrr^ = -i:i{L^I'^) , m = 0,l,2,.... (25) 
J m ^ 

where = (ri/m) Res (L™^") is called the conserved density of the charge Hm- For the 
KdV equation (PT|) we have the nontrivial charges 

r 1 

Hi = dx —u , 
7 2' 

Hi = dx -V? , 
7 8' 

= j dx^{2u^-ul) , 

I . (26) 



Another interesting feature of integrable systems is that they are bi-Hamiltonian 
i.e., these systems are Hamiltonian with at least two compatible Hamiltonian structures 
(without loss of generality we will consider from now on the KdV hierarchy instead of the 
generalized KdV hierarchy ()17j) ) 

^ = KM = {U, = Pl^^ = {U, ff42 = ^^2^ . (27) 

Otn OU OU 

Then, it can be shown that the charges are in involution with both Hamiltonian 
structures 

Hm]i = {Hn, Hm}2 = , n, m = 0, 1, 2, . . . . 
For the KdV equation we have 

6Hi 6Ho 
ut = T^i^— = V2^— , 

OU OU 

with 

V2 = d^ + 2{ud + du) . (28) 



Introducing the recursion operator 

R = V2V^' , 

it follows from ()27|) that the hierarchy of equations can be generated by 

Kn+l = ■ 

Also, using 

where i?^ is the adjoint of R, a recursion scheme for the Hamiltonians if„ can be written 
as 



6u 6u 

For a bi-Hamiltonian system of evolution equations, Ut = Kn[u], a natural Lax description 
is given by [22] 

— = [R, Bn] , 



dt 

where, we can identify 



Bn = K. 



Here K'^ represents the Frechet derivative of Kn, defined by 

-^nM V = ^ Kn[u + ev] 

It can also be shown that 

and the conserved charges are given by 

H2n+i OC Tr [r^^ , n = 0, 1, 2, . . . . (29) 
For the KdV equation (|2H) with bi-Hamiltonian structures (j28|) we have 



R = V2V^^ =d'^ + 2u + 2dud-^ , 

5 = ^9(9^ + 6m), (30) 
and the conserved charges also follows from ()29|1 . 

6 The Program 

In this section we describe the procedures in the program PSEUDO used to perform the 
calculations with PDO in integrable models. In the next section we present explicit exam- 
ples. Besides the procedures already introduced in section |21 and El to handle streams we 
have the following procedures to handle the PDO streams: 



Ifi 



Procedure name: variables 

Feature: The procedure variables defines the dependent variables, i.e., the evolution fields 
and the name of independent variables, avoiding the argument list in the input and output. 

Calling sequence: 

> variables ( [ul,u2, . .] , [x,tl,t2, . . .] ) ; 

Parameters: 

ul ,u2, . . . - name of the dependent variables or evolution fields 
x,tl ,t2 ... - name of the independent variables 

Examples: 

> read "pseudo .mpl" ; 

> variables ( [u,v] , [x,t] ) ; 



Procedure name: switch_diff 

Feature: Helps the visualization of derivatives and reduces the size of formulas. When 
activated all derivatives, for instance, of m = u{x, t) will be displayed and so 

on. This procedure is adapted from a similar one in the MAPLE package PDEtools 

Calling sequence: 

> switch_dif f (key) ; 

Parameters: 

key - the string "on" to switch the pretty derivative on and "off" to turn it off. 
Examples: 

> expr:=diff (u,x$3)/(diff (u,x)+diff (v,x$2))"2; 



u, V 




dx' 



expr := 




2 



> switch_dif f ("on 



); 



'on' 



> expr; 
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> dif f (expr ,x) ; 



'^xxxx 2 Uxxx \Uxx ~l~ 'Vxxx) 

{Ux+VxxY {Ux + VxxY 

We implement PDO, such as P in as streams with lazy evaluation. In this way, 
no truncation order needs to be given, terms are automatically produced on demand. The 
terms of the PDO P are represented as a sequence stream^, starting with the term with 
highest degree, 

P„a^ Pn-ld^-\ .... 

Procedure name: print_pseudo 
Feature: Prints the PDO stream. 
Calling sequence: 

> print_pseudo(p,n) ; 

> print_pseudo (p) ; 

> print_pseudo(p, [n] ) ; 

Parameters: 

p - a stream representing the PDO to be printed 

n - a integer corresponding to the number of terms to be displayed in the sum 

Observation: print_pseudo(p,n) prints the stream, just for visualization, as a PDO with 
n terms. If just the parameter p is passed to the procedure then the number of printed 
terms will be the value set to the global variable Order. print_pseudo(p, [n] ) prints, for 
further manipulation, the first n terms of the stream as a function. 

Procedure name: d 

Feature: Creates PDOs given by the Leibniz rules (0) and (jH)). It returns a delayed stream. 
Calling sequence: 

> d(n,f); 

> d(n); 

> d(f); 

> d([n] ,f); 

^In jB] PDO are represented as a coefficient stream. 
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Parameters: 

n - a positive or negative integer representing the order of the PDO to be created 
f - a function 

Observation: d(n,f ) stands for d"'f. d also accepts one argument. If n is absent then d(f ) 
represents the function / as a stream and if f is absent then d(n) represents the pure PDO 
d"- as a stream. d([n] ,f) stands for (c?"/), i.e., if n > it represents as a stream the 
function 

/) = fxx^^ 1 

n 

otherwise, for n < 0, it represents^ as a stream the integral 

—n 

The procedure d prints a sequence of integers as terms are forced by the first time (we will 
not show them in the output of the examples that follows) . 

Examples: 

> d(3,u); 

"Delayed" 

> print.pseudoC/o) ; 

u + 3 Ux d"^ + S Uxx d +Uxxx + --- 
We can also use display_stream to visualize a stream 

> display_streain(d(3,u) ,20) ; 



d^u, 3d^Ux, 3uxxd, Uxxx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 

However, display_stream does not order d, in fact d in any PDO output expression is 
always meant to be at the most right position. This example also shows us that we do not 
make distinction between infinite and finite PDO when we represent them as streams. 

> print_pseudo(d(-l,v) ,3) ; 



vd ^ — v^d ^ + v^,^ d ^ + . . • 

''We mean an algebraic representation and not an analytical representation. The last one is given by 
nonlocal integrals. 
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> Order: =5: 

> print _pseudo(d(-l ,v) ) ; 

vd"^ — Vxd""^ + Vxx — Vxxx d"^ + 

> print_pseudo(d(-3)) ; 

d-' + ... 

> print_pseudo(d(u''2) ) ; 



> print_pseudo(d(u"2) , 100) ; 



> display_streain(d(u''2) , 10) ; 

u^, 0, 0, 0, 0, 0, 0, 0, 0, 

> print_pseudo(d( [-2] ,u)) ; 

r I 

u dx dx + 



> II:=print_pseudo(d([-2] ,u) , [1]) ; 

> print_pseudo(d( [2] ,u)) ; 

> DD:=print_pseudo(d([2] ,u) , [1]) ; 

DD := Uxx 

> print_pseudo(d([l] ,11)) ; 

90 



> print_pseudo(d([-l] ,DD)) ; 

+ ... 

With the procedure d and the following ones we can construct any PDO. 

Procedures names: multiply_pseudo or &*, add_pseudo or &+, subtract_pseudo or &- 
Feature: Multiplies, adds or subtracts any number of PDOs returning a delayed stream. 
Calling sequence: 

> multiply_pseudo(pl,p2,p3, . . .) ; 

> &*(pl,p2); 

> pi &* p2; 

> add_pseudo(pl,p2,p3, . . . ) ; 

> &+(pl,p2); 

> pi &+ p2; 

> subtract_pseudo(pl,p2,p3, . . . ) ; 

> &-(pl,p2); 

> pi &- p2; 

Parameters: 

pl,p2,p3, . . . - the streams representing the PDOs to be multiplied, added or subtracted 
Procedure name: scalar_multiply_pseudo or &. 

Feature: Multiply a PDO by any number of functions returning a delayed stream. 
Calling sequence: 

> scalar_inultiply_pseudo(f 1 ,f 2, . . . ,p) ; 

> &.(f,p); 

> f &. p; 

Parameters: 

f l,f2, . . . ; f - the functions multiplying the PDO 
p - a stream representing the PDO 

Observation: For just two arguments the procedures &*, &+, &- and &. can be used as an 
unary prefix operator or as an infix binary operator. We must be very careful and observe 
the order of precedence of the & operators when used as binary operators, they are all left- 
associative and &* has the lowest binding strength: &+, &-, &., &*. Parentheses should be 
used to avoid dubioTis inputs. Also, there is no operator overloading, therefore, care must 
be taken in order to not mix streams with functions. 

Examples: 
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> print_pseudo((v"2 &. d(l,u))); 

Error, (in force) 'delay_record' does not evaluate to a module 

> print_pseudo( ( (v~2) &. d(l,u))); 

v'^ Ud + Ma; + ... 

> print_pseudo (-5 &. d(2)); 

Error, (in force) 'delay_record' does not evaluate to a module 

> print_pseudo((-5) &. d(2)); 

-5d^ + ... 

> print_pseudo( d(l) &* d(l) &+ d(2) &- (2 &. d(2))); 

-d^ + d^ + ... 

> print_pseudo( (d(l) &* d(l)) &+ d(2) &- (2 &. d(2))); 

zero + ... 

For the example given by equations ()12|) 

> P:=d(l) &+ ((u"2) &. d(-l)) : 

> Q:=d(2) &+ (u &. d(-l)): 

> print_pseudo(P &* Q,8); 

d^ + d + u + Ux d~^ + d^"^ — Ux d^^ + u^x d~^ + ... 

The multiply_pseudo procedure uses the strategy described in Section HI This is 
because implementing the product of PDO's as we have done for power series in Sectional 
is not straightforward. Terms of PDO acting on terms of another PDO produce another 
PDO, i.e, a infinite stream instead of a function. The other algebraic procedures are slightly 
variations of the operations implemented in Section El for power series. 

Procedure name: negate_pseudo 

Feature: Multiplies all terms of a PDO by —1. 

Calling sequence: 

> negate_pseudo(p) 

99. 



Parameters: 
p - the stream representing the PDO 

Observation: negate_pseudo(p) is a syntactic sugar for ((-1) &. p). 
Procedure name: nth_root_pseudo 

Feature: Calculates the n*^ root of the PDO P in the PDO v^. 
Calling sequence: 

> nth_root_pseudo(p,n) 

Parameters: 

p - the stream representing the PDO P 
n - a integer representing the root n 

Observation: n should divide the order of the PDO P. n can be negative. 
Examples: 

Using the operator (IHj), 

> L:=d(2) &+ d(u) : 

> A:=nth_root_pseudo(L,2) :print_pseudo(A,6) ; 

d +-ud --u^d +(^y + g"^^)" + \^uux- —Uxxx] d + .. 
which is the square root in (jl5|) . 

> B : =nth_root_pseudo (L , -2) : print_pseudo (B , 6) ; 

,-i 1 j-3 3 /3m2 7 \ /15 15 
a ud H — M^d + Uxx \ d + — Uxxx uux \ d + 



> print_pseudo(A &* A, 10); 



> print_pseudo(A &* B,10); 



d"^ + U + ... 



1 + 



Procedure name: invert\_pseudo 



Feature: Gives the inverse of a PDO. 
Calling sequence: 
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> invert_pseudo(p) 

Parameters: 
p - the stream representing the PDO 

Observation: invert_pseudo(p) is a syntactic sugar for nth_root_pseudo(p,-l). 
Procedure name: mnth_root_pseudo 

Feature: Calculates the m*^ power of the n*^ root of the PDO P in (fT!^ . the PDO {\fpy 
Calling sequence: 

> mnth_root_pseudo(p,m,n) 

Parameters: 

p - the stream representing the PDO P 
m - a integer representing the power m 
n - a integer representing the root n 

Observation: For m = 1 it behaves like nth_root_pseudo. m and/or n can be negative. 
Examples: 

> C:=mnth_root_pseudo(L,3,2) :print_pseudo(C,6) ; 

j3 3 ^ 3 fSu^ ^ \ ,-1 f ^ 1 \ 

d'' + -ud + -Ux + { —— + -Ua;x ] d + { --UU^ - —Uxxx ] d + ... 



2 4 V 8 8 / V 8 16 

which is the cubic root in (jl5|) . Since nth_root_pseudo and mnth_root_pseudo accept 
negative arguments we can use these procedures to divide PDOs: 

Procedure name: &/ 

Feature: Divides two PDOs returning a delayed stream. 
Calling sequence: 

> &/(pl,p2); 

> pi &/ p2; 

Parameters: 

pl,p2 - the streams representing the PDOs to be divided 

Observation: k/ (pi ,p2) is a syntactic sugar for (pi &* nth_root_pseudo(p2,-l)). 
Examples: 

> print_pseudo((L &/ C) &- invert_pseudo(A) , 10) ; 
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zero + ... 

We do not work with inert PDOs, that means that any PDO is expanded, for instance d~^u 
yields 

> 'd(-l,u) '=print_pseudo(d(-l,u) ,5) ; 



d{—l,u)=ud ^ — Uxd '^ + Uxxd ^ — Uxxxd ^ + Uxxxxd ^ + ... 
Exphcit formulae for the coefficients of PDO resulting from some operation, such as 

1 - ^--1/25-1^-1/2 , 



-u 



du + ud 2 

are not easy to obtain in general and are avoided in our program. However, we can still 
check the validity of this equality. For the operator OP = du + ud 

> OP:=d(l,u) &+ (u &. d(l)) : 

> OP_invert :=invert_pseudo(OP) : 

> print_pseudo(0P_invert,3) ; 



1 d 1 1 Uxd 1 {—3ux'^ + 2uUxx) d ^ 

2 M 4 u2 8 

> print_pseudo(OP_invert &* OP, 10); 



1 + 

> 0P_invert_exact:=(u~(-l/2)/2) &. d(-l ,u" (-1/2) ) : 

> print_pseudo (OP_invert_exact &* OP, 10); 



1 + 

> print_pseudo(OP_invert_exact &- OP_invert , 10) ; 



zero + ... 



Procedure name: residue_pseudo 
Feature: Gives the residue fl22j) of the PDO. 

Calling sequence: 

> residue_pseudo(p) ; 



Parameters: 

p - a stream representing the PDO 
Examples: 

> residue_pseudo(C) ; 

3^2 1 



Procedure name: project_pseudo 

Feature: Returns the part of the PDO as indicated in |TH)). The result is returned as a 
stream. 

Calling sequence: 

> project_pseudo(p,d) 

Parameters: 

p - a stream representing the PDO 

d - a integer representing the degree of the projection 

Examples: 

> E : =pro j ect_pseudo (C , 0) : pr int_pseudo (E) ; 

,33^3 
d + - ud + -u^ + ... 

Procedure name: commutator.pseudo 

Feature: Calculates the commutator of two PDO returning a stream as the resulting PDO. 
Calling sequence: 

> comniutator_pseudo(pl ,p2) 
Parameters: 

pi ,p2 - the streams representing the PDOs 
Examples: 

> F : =coinmutator_pseudo (L , E) : pr int_pseudo (F) ; 

_3 _ 1 

2 ^ 4 

Procedure name: dif f erentiate_pseudo 

Feature: Differentiates each term of a PDO with respect with some variable returning a 
stream as the resulting PDO 

Calling sequence: 

9fi 



> different iate_pseudo(p,x) 

Parameters: 

p - a stream representing the PDO to be differentiated 
X - the variable of differentiation 

Procedure name: equations_pseudo 

Feature: Prints, just for visuahzation, a number of equations resulting in taking each coef- 
ficients of a stream and making it equal to zero. Useful to print the equations produced by 
Lax equations. 

Calling sequence: 

> equations_pseudo(p,n) 

Parameters: 

p - a stream representing a PDO 

n - a integer representing the number of equations to be printed 
Examples: 

> equations_pseudo(dif f erentiate_pseudo(L,t) &+ F,l); 



Ut - -UUx - - Uxxx = 

This is of course the KdV equation ()21|) that follows from the Lax pair of operators ()20j) 
after using the Lax equations (fTB)). 



7 Examples 

In this section we present an interactive MAPLE session with the Input /Output data of 
some classical examples as well some results obtained recently with the aid of the program 
PSEUDO. We start with 

> read "pseudo .mpl" : 

> switch.diff ("on") : 

> Order: =10: 

7.1 The KdV Hierarchy 

The classical KdV hierarchy of equations follows from the Lax operator ()14|). 

> variables ( [u] , [x,t] ) : 

> L_kdv:=d(2) &+ d(u) ; 
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L_kdv := "Delayed" 

> print.pseudo (L_kdv) ; 

+ U + ... 

The square root (fT5|l can be calculated from 

> print_pseudo (mnth_root_pseudo (L_kdv ,1,2)); 



1 ,_i 1 ._o / ^2 1 \ /3 1 

a 

7 11 o 1 



d H — ud ^ Uxd ^+( 1 — Uxx I d ^ + { - uux Ux 

2 4 V 8 8 / V8 16 



~l~ ( -|^g ^^xx 22 ~^ 22 ^^^^^ J ^ 

15 2 15 15 _1_ \ 6 

I I Hi Hi ^ I W W XXX I X 3!)IC < xxxxx I 

32 32 16 64 / 

/85 2 55 2 31 37 1 91 2 5m^\ 7 

+ 1 U Ua; H -U Uxx UU 

XXXX irnr\ ^XXX ~l~ ^XXXXXX ^-1/-, ^XX I '-^ 

.64 64 64 32 128 128 128 

175 175 35 3 175 2 1 175 3 

~l~ I 22 ^ ~l~ "^xxxx ~l~ W Y28 ^ ^^'^^ 256 T28 

245 63 

I XX XXX I xxxxx I I • • • 

128 128 



The Lax equation 



dL 



KdV _ r/ rn/2 



[(^K^dv) + 5 -t^KdVj ! 
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will produce nontrivial equations for n odd. The first equations are 

> for n from 1 to 8 by 2 do equations_pseudo(dif f erentiate_pseudo(L_kdv,t) 
&+ commutator_pseudo (L_kdv , pro j ect_pseudo (mnth_root_pseudo (L_kdv , n , 2) , 0) ) 
,10) end do; 

Ut — Ux = 





3 




1 


= 




Ut 

2 


U Ux 


^ ^xxx 


15 p 


5 


5 




1 


— 

8 


HL ^t_Jj ^t_Jj XXX 

8 


~ 4 


'^xx 


— u 
16 ■ 




35 


7 




35 


'^xxx ~ 


^i-ll '\JL rj* HL XX 

8 


~ 32 


xxxxx 


"32' 



35 3 35 2 35 7 35 3 35 21 



16 32 8 32 32 32 32 

1 

„ . '^xxxxxxx ~l~ 

64 

The first nontrivial odd conserved densities are 



98 



> for n from 1 to 8 by 2 do residue_pseudo(mnth_root_pseudo(L_kdv,n,2)) 
end do; 



u 
2 

3^^ 1 
8 8 

5u^ 5 5 2 1 

16 16 32 32 

35 2 7 35 2 7 35 21 2 1 

~ Uj Uj 2; | U Uj XXX I T XX I T XXXX I | "Z^ | xxxxxx 

64 32 64 64 128 128 128 

Assuming that u ^ sufficiently fast as |x| ^ oo we obtain from ()25|) (after integration by 
parts and up to normalization constants) the conserved charges ()26p. 

Ahernatively, we can obtain ^3] (up to normahzation constants) all results for the 
KdV hierarchy using the recursion operator (j30|) 

> R_kdv:=d(2) &+ d(4*u) &+ ((2*diff (u,x)) &. d(-l)); 

R_kdv := "Delayed" 

> print_pseudo(R_kdv) ; 

d"^ + 4:U + 2uxd'^ + ... 

> print_pseudo (mnth_root_pseudo (R_kdv ,1,2)); 



d +2ud ^-2u^d '■^ + Auuxd'^ + {Au^ -2ux^ - 4:UUxx)d~^ 

+ ( — 24m^ Ux + AUxUxx + 4:UUxxx) d~^ 

+ (52 u Ux^ + 40 Uxx — 2 Uxx^ — 10 — 4 m Uxxxx — 4 Uxxx) d~^ 

I ^ 4 ?X X U^ xxxx 3 6 Uj X 192 Uj Uj ^ U^ XX I 4 U^ U^ xxxxx I 4 xx xxx I 120 "XIj U ^ 6 Uj Uj xxx ) 

> for n from 1 to 8 by 2 do 

equations_pseudo (dif f erentiate_pseudo (R_kdv,t) &+ 

commutator_pseudo (R_kdv , pro j ect_pseudo (mnth_root_pseudo (R_kdv , n , 2) , 0) ) 
,10) end do; 
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2 Uxx 



= 
. = 



4:Ut-24:UUx- 4: U^XX = 

2 Utx 12 Ux 12 U Uxx 2 Uxxxx — 



Aut-4: Uxxxxx - 120 U Ux-4QU Uxxx - 80 Ux Uxx = 
Ifr — 120 UUt — 60 Wttt — 60 li-rr " 



4:Ut — 560 lia; — 280 u^ Uxxx — 1120 u Ux Uxx — 56 w Uxxxxx — 280 li^;^ — 168 Ux Uxxxx 

—840 It Ux — SAO UUx Uxxx — ^^OUx Uxx — ^^'^UxUxxxxx — '^SOu Uxx — ^^OuUxx 

> for n from 1 to 8 by 2 do residue_pseudo(iiinth_root_pseudo(R_kdv,n,2)) 
end do; 

2u 

6 + 2 Uxx 

20 + 20 u u^^ + 10 vj + 2 u^xx 

140 U Ux + 56 Ma; Uxxx + 140 + 28 M Uxxxx + 70 ti** + 42 Majj;^ + 2 tia;a;a;a:a;a; 

7.2 The KP Hierarchy 

In order to treat the KP hierarchy we need to construct a generic PDO, Lkp- This is easily 
done by using the mcLke_streain procedure 

> variables ( [U] , [x , t [1] , t [2] , t [3] ] ) : 

> L_KP:=add_pseudo(d(l) ,maJje_streain(proc(k) 

option remember; 
U[k] *d^(-k); 
end proCjD) : 

> print_pseudo(L_KP) ; 



The Lax equation 



dL 



KP 



dtr 



KP)+ > -^KPJ , 



n= 1,2,... 



will produce the equations in the hierarchy. For n = 1, 2, 3 we obtain 

> equations_pseudo(differentiate_pseudo(L_KP,t [1] ) &+ 
commutator_pseudo (L_KP , pro j ect_pseudo (L_KP , 0) ) , 3) ; 





-Us. 


= 






= 






= 



> equations_pseudo(dif f erentiate_pseudo(L_KP,t [2] ) &+ 
commutator_pseudo(L_KP,project_pseudo(L_KP &* L_KP,0)),3); 

U2t, - 2 f/i f/;, - 2 U3, - U2.., = 
U lt2 — '^U 2x — U ixx =0 

Ust, -AU2U1X + 2 Ui Uixx -2U^x- Usxx = 

> equations_pseudo(diff erentiate_pseudo(L_KP,t [3] ) &+ 
commutator_pseudo(L_KP,project_pseudo(L_KP &* L_KP &* L_KP,0)),3); 

Ulta — QUiUix — SUsx — 3U2XX — U ixxx = 
U2t3 — QU2U1X ~ QU1U2X — 3U^x — SUsxx — U 2xxx = 
^3*3 -~ QU^U^x ~ ^U2U 2x + ^U2U pxx + ^UiU 2xx —^Usx ~^UiUsx ~^U^xx ~U 3xxx = 



7.3 The Harry Dym Hierarchy 

The Harry Dym equation 

1 15 Wx^ — ISWxW Wxx + 4:W^Wx 



Wt = (W ^l'^)xxx 



w 



(7/2) 



and its hierarchy can be obtained from the Lax representation 

1 



dt 



w 



-2" 



2n + l 

-^HD I ^2 ' -^HD 



, n = 0,l,2,..., 



Therefore, the first equations are 



> assume(w(x,t)>0) : 

> variables ( [w] , [x , t] ) : 

> L_hd:=(l/w) &. d(2) :print_pseudo(L_hd) ; 



h ... 

w 



> for n from to 3 do equations_pseudo(dif f erentiate_pseudo(L_hd,t) &+ 
commutator_pseudo (pro j ect_pseudo (iiinth_root_pseudo (L_hd , 2*n+l , 2) , 2) , 
(2"n) &. L_hd),5) end do; 



2 

1 SWtW^ / + 15 Wx^ — IS WxWWxx + 4: w'^ Wx^^ ^ 



w 



1 



128 



4096 



128 w V ^ / + 3465 - 6930 w Wx^ Wxx + 2520 w'^ Wx Wxx^ + 1820 Wx 

17' 

560 Wxx Wxxx - 320 Wx Wxxxx + 32 W^xixa: 



■19\ 

4096^/;twV ^ / +360192 Wx Wxx Wxxxx ~ 3052896 W^ Wx^ Wxxx Wx 




+ 225792 w^ Wx w^xx - 10720710 w Wx^ Wxx + 2942940 w^ wj^ Wxxx 
— 591360 w^ Wx Wxxxx — 32256 w^ Wxxx Wxxxx — 21504 w^ Wxx Wxxxxx 
+ 310464 Wx^ Wxxx + 87360 iw^ Wx Wxxxxx - 8960 Wx w^ Wxxxxxx 
+8072064 w'^ Wx^ Wxx^ - 1397088 w^ Wxx^ Wx + 512 w^ Wxxxxxxx 

\ I i-) 

+ 3828825 Wx'' / w V ^ / = 
The charges follows from 



H-in+i) = TrLjjD , n = 0, 1, 2, . . . , 
and the first conserved densities are 



> for n from to 3 do residue_pseudo(innth_root_pseudo(L_hd,2*n-l,2)) 
end do; 



w 

1 5 Wx'' — 4 w w^x 
32 Tb 

1 —1155 Wx'^ + 1848 Wx'^ w Wxx — 448 tv"^ Wxxx — 336 Wxx^ + 64 Wx 



2048 /II 
1 



(42944 w^^^ - 425425 Wx^ - 569712 Wx^ Wx^ + 181632 Wxx Wxxx Wx 



65536 

+ 52800 Wx'^ Wxxxx ~ 6912 Wx + 1021020 Wx'^ w Wxx — 8832 w'^ Wxxx^ 

ir 



274560 w^^^ + 512 14592 tu"^ Wxzxz) 



2 



After integration by parts and up to normalization constants the conserved charges can be 
written as 



H-2 = jdx\w-''Wx, 



= I dx-^r^w-^'/^wt-w-y^wi 



4 



6 Vl6 



The Hunter-Zheng equation 

also belongs to the Harry Dym hierarchy. Its Lax representation is given by 

dLuD 



dt 

where 



= -2[S,Lhd], 



B = ^(d-^w)d+^d-\d-^w)d^ . 



In fact, we can check this result 



> B:=((l/4) &. (d([-2],w) &* d(l))) &+ ((1/4) &. (d(-l) &* d([-2],w) 
&* d (2) ) ) : print _pseudo (B , 4) ; 



Iff If 1 1 

-/ wdxdxd — - w dx + -w d ~^ — -Wx d ~^ + 

> equations_pseudo(dif f erentiate_pseudo(L_hd,t) &+ 

commutator_pseudo(scalar_multiply_pseudo(2,B) ,L_hd) , 1) ; 



'Wt + 2 J wdxw + J J w dx dx 



8 Conclusion 

Using lazy evaluation and streams we have described our set of routines to perform calcu- 
lations with integrable models using PDOs. In order to make the techniques used available 
for a wide audience we have introduced in detail the procedures to incorpore lazy evaluation 
as well streams in MAPLE. 

Our program provides just a basic set of operations in a way that more complicated 
procedures can be easily constructed. Our program just works for one dimensional PDOs 
. We intend to generahze this program to incorporate PDO with more general coefficients, 
such as matrix. A version to incorporate supersymmetric PDOs, among others generaliza- 
tions, is under development. 
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